function N = reference_basis_function_2d(basis_type)

switch basis_type
    case "P1"
        N = cell(3,1);
        N{1} = @(xi,eta) 1 - xi - eta;
        N{2} = @(xi,eta) xi;
        N{3} = @(xi,eta) eta;

    case "P1b"
        N = cell(4,1);
        N{1} = @(xi,eta) 1 - xi - eta;
        N{2} = @(xi,eta) xi;
        N{3} = @(xi,eta) eta;
        N{4} = @(xi,eta) 27*N{1}(xi,eta).*N{2}(xi,eta).*N{3}(xi,eta);

    case "P2"
        N = cell(6,1);
        N{1} = @(xi,eta) 1 - 3*xi + 2*xi.^2 - 3*eta + 2*eta.^2 + 4*xi.*eta;
        N{4} = @(xi,eta) 4*xi - 4*xi.^2 - 4*xi.*eta;
        N{2} = @(xi,eta) -xi + 2*xi.^2;
        N{5} = @(xi,eta) 4*xi.*eta;
        N{3} = @(xi,eta) -eta + 2*eta.^2;
        N{6} = @(xi,eta) 4*eta - 4*eta.^2 - 4*xi.*eta;

    case "Q1"
        N = cell(4,1);
        N{1} = @(xi,eta) 0.25*(1 - xi - eta + xi.*eta);
        N{2} = @(xi,eta) 0.25*(1 + xi - eta - xi.*eta);
        N{3} = @(xi,eta) 0.25*(1 + xi + eta + xi.*eta);
        N{4} = @(xi,eta) 0.25*(1 - xi + eta - xi.*eta);

    case "Q1b"
        N = cell(5,1);
        N{1} = @(xi,eta) 0.25*(1 - xi - eta + xi.*eta);
        N{2} = @(xi,eta) 0.25*(1 + xi - eta - xi.*eta);
        N{3} = @(xi,eta) 0.25*(1 + xi + eta + xi.*eta);
        N{4} = @(xi,eta) 0.25*(1 - xi + eta - xi.*eta);
        N{5} = @(xi,eta) 256.*N{1}(xi,eta).*N{2}(xi,eta).*N{3}(xi,eta).*N{4}(xi,eta);

    case "Q2"
        N = cell(9,1);
        N{1} = @(xi,eta) 0.25*(xi.*eta - xi.*eta.^2 - xi.^2.*eta + xi.^2.*eta.^2);
        N{5} = @(xi,eta) 0.5*(-eta + eta.^2 + xi.^2.*eta - xi.^2.*eta.^2);
        N{2} = @(xi,eta) 0.25*(-xi.*eta + xi.*eta.^2 - xi.^2.*eta + xi.^2.*eta.^2);
        N{6} = @(xi,eta) 0.5*(xi + xi.^2 - xi.*eta.^2 - xi.^2.*eta.^2);
        N{3} = @(xi,eta) 0.25*(xi.*eta + xi.*eta.^2 + xi.^2.*eta + xi.^2.*eta.^2);
        N{7} = @(xi,eta) 0.5*(eta + eta.^2 - xi.^2.*eta - xi.^2.*eta.^2);
        N{4} = @(xi,eta) 0.25*(-xi.*eta - xi.*eta.^2 + xi.^2.*eta + xi.^2.*eta.^2);
        N{8} = @(xi,eta) 0.5*(-xi + xi.^2 + xi.*eta.^2 - xi.^2.*eta.^2);
        N{9} = @(xi,eta) 1 - xi.^2 - eta.^2 + xi.^2.*eta.^2;

    otherwise
        error("Invalid basis type.");
end

end